%% 
clearvars;
%% read the data
root=('D:\Documents2\Matlab\data_analysis\Sn\171219 MPI mK meas\critical field 180808\'); %
filenames={'2Sn_15PbTe';'3Sn_6PbTe';'3Sn_8PbTe';'3Sn_12PbTe'};
ii=4;
 file1=[root,filenames{ii},'_op.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f','commentStyle','T');
                            % Tmc(K)  Hc2(T) 
        data1=data1';     
        fclose(fid);
        
  file1=[root,filenames{ii},'_ip.dat'];
        fid=fopen(file1);
        data2 = textscan(fid,'%f %f','commentStyle','T');
                            % Tmc(K)  Hc2(T) 
        data2=data2';     
        fclose(fid);
        
       
figure   
plot(data2{1},data2{2},'ko');
 hold on
%% Fit the in-plane data 2D Ginzburg-Landau
s2DGL = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[0.1],...
               'Upper',[20],...
               'Startpoint',[5]);
f2DGL = fittype('a*x','options',s2DGL);
index=find(data2{1}>=0.4*max(data2{1}));

tt=data2{1}/max(data2{1});
[fitcoef,gof1]=fit(sqrt(1-tt(index)),data2{2}(index),f2DGL);
c_2DGL=coeffvalues(fitcoef);
T1=linspace(0,1,201)';
Hc1=c_2DGL*sqrt(1-T1);
plot(T1*max(data2{1}),c_2DGL*sqrt(1-T1),'-r');
 hold on
 
        
%% Fit the in-plane data with KLB
% supplement of DOI: 10.1038/NPHYS3580, Saito, Ising SC in MoS2
sKLB = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[0],...
               'Upper',[10],...
               'Startpoint',[0.001]);
fKLB = fittype('-psi(0,0.5+A*x)+psi(0,0.5)','options',sKLB);

X=data2{2}.^2./data2{1};
Y=log(data2{1}/max(data2{1}));
index=find(data2{1}>=0.4*max(data2{1}));

[fitcoef,gof2]=fit(X(index),Y(index),fKLB); 
c_KLB=coeffvalues(fitcoef);
  Hc_2_T=(0:0.01:1000)';
  lnT=-psi(0,0.5+c_KLB*Hc_2_T)+psi(0,0.5);
  T2=exp(lnT)*max(data2{1});
  Hc2=sqrt(T2.*Hc_2_T);
plot(T2,Hc2,'-g');
 if ii<4
    axis([0 0.55 -0.2 4]);
 else
     axis([0 1.2 -0.2 4.3]);
 end
hold on
     

%% Fit out-of-plane critical field with Gurevich
figure
plot(data1{1},data1{2},'ko');

if ii<4
    axis([0 0.55 -0.01 0.2]);
 else
     axis([0 1.2 -0.03 0.6]);
 end
hold on


if ii==4
    T0=0.01;
else
    T0=0.01;
end
index=find(data1{1}>=T0);
lnt2=log(data1{1}/max(data1{1}));
h2=data1{2}./data1{1};    

options.MaxIter=300;
options.RobustWgtFun='welsch';
p_Gurevich_op=nlinfit(h2(index),lnt2(index),@Gurevich07,[2 0.4 0.7 0.5 0.01],options);

Hc_T=(0:0.005:1000)';
lnT=Gurevich07(p_Gurevich_op,Hc_T);
T3=exp(lnT)*max(data1{1});
Hc3=T3.*Hc_T;

plot(T3,Hc3,'-k');




%% End




